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Abstract 

We consider the Gross-Pitaevskii (GP) equation in the presence of pe- 
^vj ■ riodic and quasiperiodic superlattices to study cigar-shaped Bose-Einstein 

^ ' condensates (BECs) in such potentials. We examine spatially extended 

fT^ , wavefunctions in the form of modulated amplitude waves (MAWs) . With 

\^ ■ a coherent structure ansatz, we derive amplitude equations describing the 

("^ ' evolution of spatially modulated states of the BEG. We then apply second- 

^^D , order multiple scale perturbation theory to study harmonic resonances 

^— >^ ' with respect to a single lattice wavenumber as well as ultrasubharmonic 

'^ ' resonances that result from interactions of both wavenumbers of the su- 

perlattice. In each case, we determine the resulting system's equilibria, 
which represent spatially periodic solutions, and subsequently examine 
the stability of the corresponding solutions by direct simulations of the 
P^ . GP equation, identifying them as typically stable solutions of the model. 

We then study subharmonic resonances using Hamiltonian perturbation 
theory, tracing robust, spatio-temporally periodic patterns. 
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1 Introduction 

At very low temperatures, particles in a dilute bose gas can occupy the same 
quantum (ground) state, forming a Bose-Einstein condensate (BEC) 0H1 I^H 



I3(J[ I17| . which appears as a sharp peak (over a broader distribution) in both 
coordinate and momentum space. As the gas is cooled, condensation (of a large 
fraction of the atoms in the gas) occurs via a quantum phase transition, emerging 
when the wavelengths of individual atoms overlap and behave identically. Atoms 
of mass m and temperature T constitute quantum wavepackets whose spatial 
extent is given by the de Broglie wavelength 



^db = \ — p^ , (1) 

which represents the uncertainty in position associated with the momentum 
distribution PUJ (where h is Planck's constant and feg is Boltzmann's constant). 
The atomic wavepackets overlap once atoms are cooled sufficiently so that Xdt is 
comparable to the separation between atoms, as bosonic atoms then undergo a 
quantum phase transition to form a BEC (a coherent cloud of atoms). Although 
a condensate constitutes a quantum phenomenon, such "matter waves" can often 
be observed macroscopically, with the number of condensed atoms N ranging 
from several thousand (or less) to several million (or more) |21) . 

BECs were first observed experimentally in 1995 in dilute alkali gases such 
as vapors of rubidium and sodium 01 122| . In these experiments, atoms were 
confined in magnetic traps, evaporatively cooled to a fraction of a microkelvin, 
left to expand by switching off the confining trap, and subsequently imaged 
with optical methods. A sharp peak in the velocity distribution was observed 
below a critical temperature, indicating that condensation had occured [as the 
alkali atoms were now condensed in the same (ground) state]. Under the typi- 
cal confining conditions of experimental settings, BECs are inhomogeneous, so 
condensates arise as a narrow peak not only in momentum space but also in 
coordinate space. 

The macroscopic observability of the condensates in coordinate and momen- 
tum space has led to novel methods of investigating quantities such as energy 
and density distributions, interference phenomena, the frequencies of collective 
excitations, the temperature dependence of BECs, among others |2] (for com- 
prehensive reviews, the interested reader should consult Refs. 0H1EZI)- Another 
consequence of this inhomogeneity is that the effects of two-body interactions 
are greatly enhanced, despite the fact that bose gases are extremely dilute (with 
the average distance between atoms typically more than ten times the range of 
interatomic forces). For example, these interactions reduce the condensate's 
central density and enlarge the size of the condensate cloud, which becomes 
macroscopic and can be measured directly with optical imaging methods. 

BECs have two characteristic length scales. The condensate density varies on 
the scale of the harmonic oscillator length Uho — y'njirriujho) [which is typically 
on the order of a few microns], where Uho = {^x^y^zY is the geometric 
mean of the trapping frequencies. The "coherence length" (or "healing length" ) , 
determined by balancing the quantum pressure and the condensate's interaction 
energy, is x = l/-\/87r|a|Ti [and is also typically of the order of a few microns], 
where n is the mean particle density and a, the (two-body) s-wave scattering 



length, is determined by the atomic species of the condensate. Interactions 
between atoms are repulsive when a > and attractive when a < 0. For a 
dilute ideal gas, a « 0. The length scales in BECs should be contrasted with 
those in systems like superfluid helium, in which the effects of inhomogeneity 
occur on a microscopic scale fixed by the interatomic distance I21j . 

If considering only two-body, mean-field interactions, a dilute Bose-Einstein 
gas near zero temperature can be modeled using a cubic nonlinear Schrodinger 
equation (NLS) with an external potential, which is also known as the Gross- 
Pitaevskii (GP) equation. This is written [21] 

^h%^(^~^+g,\^\' + V{r}^^, (2) 

where ^ = ^(r, t) is the condensate wave function normalized to the number 
of atoms, V(r) is the external potential, and the effective interaction constant 
is go = [47r/i^a/TO][l + 0{(^)], where C = VT^I^H^ is the dilute-gas parameter 
[21ii35n7 . 

BECs are modeled in the quasi-one-dimensional (quasi-lD) regime when 
the transverse dimensions of the condensate are on the order of its healing 
length and its longitudinal dimension is much larger than its transverse ones 
|18l 1141 [T^ I21j . In this regime, one employs the ID limit of a 3D mean-field 
theory (generated by averaging in the transverse plane) rather than a true ID 
mean-field theory, which would be appropriate were the transverse dimension 
on the order of the atomic interaction length or the atomic size ^J USJ |S] . The 
resulting ID equation is |5flll21j 

ihut — —[fi^ / {2m)]uxx + glu] u + V{x)u, (3) 

where u, g, and V are, respectively, the rescaled ID wave function ("order pa- 
rameter"), interaction constant, and external trapping potential. The quantity 
|wp gives the atomic number density. The self- interaction parameter g is tunable 
(even its sign), because the scattering length a can be adjusted using magnetic 
fields in the vicinity of a Feshbach resonance |24l I34| . The manipulation of 
Feshbach resonances has become one of the most active areas in the study of 
ultracold atoms, as (for example) numerous research groups are investigating 
the intermediate regime between molecular condensates and degenerate Fermi 
gases (the so-called "BEC-BCS" crossover regime). Theoretical algorithms for 
manipulating a, such as alternating it periodically between positive and nega- 
tive values, have been developed by analogy with "dispersion management" in 
nonlinear optics. 

In forming a BEG, the atoms are trapped using a confining magnetic or op- 
tical potential V{x), which is then turned off so that the gas can expand and 
be imaged. In early experiments, only parabolic ("harmonic") potentials were 
employed, but a wide variety of potentials can now be constructed experimen- 
tally. In addition to harmonic traps, these include double-well traps (see, e.g., 
|3] and references therein), periodic lattices (see, e.g., JJj for a review), and 



super lattices 07||^ (which can be either periodic or quasiperiodic), and super- 
positions of lattices or superlattices with harmonic traps. Optical lattices and 
superlattices are created using counter-propagating laser beams, and higher- 
dimensional versions of many of the aforementioned potentials have also been 
achieved experimentally. 

The existence of quasi-lD ("cigar-shaped") BECs motivates the study of 
lower dimensional models such as Eq. Q . The case of periodic and quasiperi- 
odic potentials without a confining trap along the longitudinal dimension of the 
lattice is of particular theoretical and experimental interest. Such potentials 
have been used, for example, to study Josephson effects |3], squeezed states 03], 
Landau-Zener tunneling and Bloch oscillations |321) and the transition between 
superfluidity and Mott insulation at both the classical [SBJ E] and quantum 
pS) levels. Moreover, with each lattice site occupied by one alkali atom in its 
ground state, BECs in optical lattices show promise as a register in a quantum 
computer [SIEHI- 

In experiments, a weak harmonic trap is typically used on top of the optical 
lattice (OL) or optical superlattice (OSL) to prevent the particles from escaping. 
The lattice is also generally turned on after the trap. If one wishes to include 
the trap in theoretical analyses, then V{x) is modeled by 

V{x) = Vi cos(kix) -I- V2 cos(K2a;) -I- T^x^ , (4) 

where Ki is the primary lattice wavenumber, K2 > ki is the secondary lattice 
wavenumber, Vi and V2 are the associated lattice amplitudes, and Vh represents 
the magnitude of the harmonic trap. Note that Vi, V2, Vh, ki, and K2 can all 
be tuned experimentally, so that the external potential's length scales are easily 
manipulated. The sinusoidal terms in Q) dominate for small x, but the harmonic 
trap otherwise becomes quickly dominant. When Vh ^ V'l , V2, the potential 
is dominated by its periodic (or quasiperiodic) contributions for many periods 
[ 181 150] . BECs in OLs with up to 200 wells have been created experimentally 

In this work, we let Vh = and focus on OL and OSL potentials. Spatially 
periodic potentials have been employed in experimental studies of BECs ,221 El 
133 gamins and have also been studied theoretically [El [TUl HOI EH El Ei EOl 
031 ^5" 56 , '38', 1221 ; see also the recent reviews |22E|- In experiments reported 
in 2003, BECs were loaded into OSLs with K2 = 3ki 02] • However, there has 
thus far been very little theoretical research on BECs in superlattice potentials 
|54l 1231 rVfl I25| . In this work, we consider both periodic (rational K2/K1) and 
quasiperiodic (irrational K2/K1) OSLs. 

We focus here on spatially extended solutions rather than on localized waves 
(solitons). For BECs loaded into OSLs, the interest in such extended wavefunc- 
tions is twofold. First, BECs were successfully loaded into OSL potentials in 
recent experiments @7| (in which extended solutions were observed). Second, 
MAWs in BECs in OSLs can be used to study period-multiplied states and 
generalizations thereof j3^ jSUl O . 

On the first front, ^^Rb atoms were loaded into an OSL by the sequential 
creation of two lattice structures. The atoms were initially loaded into every 



third site of an OL. A second periodic structure was subsequently added so that 
the atoms could be transferred from long-period lattice sites to corresponding 
short-period lattice sites in a patterned loading. 

On the second front, Machholm et al. JF studied period-doubled states 
(in |mP), interpreting them as soliton trains to attempt to explain experimen- 
tal studies by Cataliottiei al. ^\, who observed superfluid current disruption 
in chains of weakly coupled BECs in OL potentials. More recently, experi- 
mental observations of period doubled wavefunctions in BECs in OL potentials 
have now been reported |^. From a dynamical systems perspective, period- 
multiplied states arise at the center of KAM islands in phase space; the location 
and size of these islands have been estimated using Hamiltonian perturbation 
theory and multiple scale analysis J3H1 HOI 1^ • 

In this study, we investigate spatially extended solutions of BECs in peri- 
odic and quasiperiodic OSLs. We apply a coherent structure ansatz to Eq. Q, 
yielding a parametrically forced DufRng equation describing the spatial evolu- 
tion of the field. We employ second-order multiple scale perturbation theory 
to study its periodic orbits (called "modulated amplitude waves" and denoted 
MAWs), and illustrate their dynamical stability with numerical simulations of 
the GP equation. We consider harmonic (1 : 1) resonances and two types of 
ultrasubharmonic resonances — resulting from, respectively, "additive" (2:1-1-1) 
and "subtractive" (2:1 — 1) interactions — all of which arise at the 0{e^) level 
of analysis. Because ultrasubharmonic resonances result from the interaction of 
multiple superlattice wavenumbers, they cannot occur in BECs loaded into reg- 
ular OLs. We then explore subharmonic resonances using Hamiltonian pertur- 
bation theory, identifying various relevant patterns including quasi-stationary 
ones (with weak amplitude oscillations) and spatio-temporally breathing ones 
(see the details below). 

We structure the rest of our presentation as follows: We first introduce 
modulated amplitude waves and use multiple scale perturbation theory to de- 
rive "slow flow" dynamical equations that describe the resonance phenomena 
under consideration. We analyze these equations and corroborate our results 
and test the stability of the MAWs with direct numerical simulations of the GP 
equation. We then examine subharmonic resonances using Hamiltonian pertur- 
bation theory and additional numerics. Finally, we summarize our findings and 
present our conclusions. 

2 Modulated Amplitude Waves 

To study MAWs, we employ the ansatz 

u{x, t) = R{x) exp [i [e{x) - ^it]) . (5) 

When these (temporally periodic) coherent structures (jSJ are also spatially pe- 
riodic, they are called modulated amplitude waves (MAWs) |16II15| . The orbital 
stability of MAWs for the cubic NLS with elliptic potentials has been studied 
by Bronski et al \1'A\ 1121 IT^ . To obtain stability information about sinusoidal 



potentials, one takes the limit as the elliptic modulus k approaches zero |36| . 
When V{x) is periodic, the resulting MAWs generalize the Bloch modes that 
occur in the theory of hnear systems with periodic potentials [5511^ 1551 ITUl 1^ . 
In this work, we extend recent studies |49[ ISU] of the dynamical behavior of 
MAWs of BECs in lattice potentials to superlattice potentials. 

Inserting Eq. (jS} into Eq. Q, equating the real and imaginary compo- 
nents of the resulting equation, and defining S := R' yields the following two- 
dimensional system of nonlinear ordinary differential equations: 

R ^ S , 
The parameter c is given by the relation 



^'(^) = 7^- (6) 

which indicates conservation of "angular momentum" ^^l- Constant phase 
solutions (i.e., standing waves), which constitute an important special case, 
satisfy c — 0. In the rest of the paper, we restrict ourselves to this class of 
solutions, so that 

R' ^S, 

S'^-^ + ^i^' + l^l-Wfl^ (T) 

We consider the case with Vh — (which implies, in practice, that the 
harmonic trap is negligible with respect to the OSL potential for the domain of 
interest) and define 

-: 2mri _ 2mg ~, , 2m^ , , .„. 

<5:=-^, £«:=-- p^, v{x):^-^V{x), (8) 

where 

V{x) = e\Vi cos(Kia;) -I- V2 cos(K2a;)] , (9) 

the parameters 5, a, and Vj are 0(1) quantities, and the lattice wavenumbers Kj 
can either be commensurate (rational multiples of each other) or incommensu- 
rate, so that the OSL can be, respectively, either periodic or quasiperiodic. We 
let K2 > 1^1 without loss of generality, so that ki is the primary lattice wavenum- 
ber. In our numerical simulations, we focus on the case K2 = 3ki which has 
been achieved experimentally |47j . 

For notational convenience, we drop the tildes from (5, q, and Vj, so that Eq. 
(Q is written 

R" + 5R + eaR^ + eR[Vi cos(Kia;) + V2 cos{k2x)] = . (10) 

In this paper, we consider the case 6 > corresponding to a positive chemical 
potential. 



3 Multiple Scale Perturbation Theory and Spa- 
tial Resonances 

To employ multiple scale perturbation theory |^ \^ , we define "slow space" 
rj := ex and "stretched space" 

^ -.^ bx ^ [I + ebi + e'^b2 + 0{e^)]x . (11) 

We then expand the wavefunction amplitude R in a power series, 

R^ Ro + eRi+e^R2 + 0{e^), (12) 

and stretch the spatial dependence in the OSL potential, which is then written 

V{0 = Vi cos{kiO + V2 cos(«;20 • (13) 

Inserting these expansions, Eq. (|10|l becomes 



[1 + b^e + b^e^ + 0{e^)y 



r/^2 



d^Ro d^Ri 



ae 



,d^R2 
de 



2e[l + 5i£ + &2e^ + 0(e3)] 



de 

d^Ro d^Ri 



+ e- 



d^Ro d^Ri 



>d^R2 



+ 0(e3) 



2d^R2 



+ 0(e3) 



Qjj2 Qj^2 Qj^2 ^ 

+ S[Ro + eRi + e^R2 + 0{e^)\ + ea [i?o + eRi + £^i?2 + ^(e^)] ^ 

+ e [i?o + eRi + e^R2 + 0(e^)] [Vi cos{kiO + V2 cos(k20] = . (14) 

To perform multiple scale analysis, we equate the coefficients of terms of 
different order (in e) in turn. At 0(1) = 0(£"), we obtain 

r)2fi>„ 

SRn = 0, 



de 



which has the solution 

i?o(C, v) = A{v) cos{VSO + Bi-q) sin(\/^e) : 



(15) 



for slowly- varying amplitudes A{ri), B{ri), equations of motion for which arise 
at 0(e). 



Equating coefficients at 0{e) yields 






SRi 



2biSA-2VdB' --aAiA^ +B^) cos{V60 
2bi5B + 2V5A' - ^aBiA^ + B^) sin(%/^0 
^[-A^ + 3B^] cos(3\/^0 + ^[-3^2 + B^] sm{3VS^) 

^— cos([ki - V^]0 + -^ cos([ki + ^S\i) 

ViB . 

sii 

2 

V2A 



([ki - \/^]e) + -^ sin([Ki + \/^]e) 



VoB 



C0s([k2 - \/^]C) + 

sin([K2 - V^]C) + 



F2A 



VoB 



■ sm 



5(N + Vs]0 

{[^2 + v^m ■ 



(16) 



For Ri{£,, if) to be bounded, the coefficients of the secular terms in Eq. H16|l 
must vanish |5H[I??]. The harmonics cos(\M^) and sin(v^O ^''6 always secular, 
whereas cos(3\/5^) and sin(3-\/^0 ^^'^ never secular. The other harmonics are 
secular only in the case of 2 : 1 subharmonic resonances |3H1 ISO]) which can 
occur with respect to either the primary (ki — 2^/5) or secondary (k2 = 2\f5) 
wavenumber of the lattice. We will consider the situation in which l|16|l is non- 
resonant and turn our attention to other resonant situations at 0(£^) that arise 
from interactions between both lattice wavenumbers. [Our 0(£^) analysis below 
can be repeated in the presence of 2: 1 resonances.] At 0(e), one obtains either 
no resonance, a long-wavelength subharmonic resonance, or a short-wavelength 
subharmonic resonance. 

Equating the coefficients of the secular terms to zero in Eq. H16|l yields the 
following equations of motion describing the slow dynamics: 



A' = -biVSB 



-^B{A^+B'), 

A{A^ 



sVs 



B'). 



(17) 



We convert p7|l to polar coordinates with A{ri) = Ccos[iy9(77)] and B{r]) = 
C sin[ip[ri)] and see immediately that each circle of constant C is invariant. The 
dynamics on each circle is given by 



viv) = "^(0) + 



b^VS 



3a 

sTs 



c 



■q. 



(18) 



We examine the special circle of equilibria, corresponding to periodic orbits of 
(|2J), which satisfies 



C^ = A-" + B^ = 



ia 



(19) 



We are interested in the O(e^) effects, which we now analyze. At this second 
order of perturbation theory, BECs in OSL potentials exhibit dynamical behav- 
ior that cannot occur in BECs in simpler OL potentials (where, for example, 
solutions of type of Eq. (|19|l straightforwardly arise \bl\). 

Equating coefficients at O(e^) yields 

-Q^ + 5R.^ -{b, + 2b,)^^ -^^-^^^^^r,- '"''"^^ - 2'^^^ - ^d^ 
- RiVi cos(kiO - R2V2 cos(k20 , (20) 

where one inserts the expressions for Rq, i?i, and their derivatives into the 
right-hand-side of ((201). 

To find the secular terms in Eq. (|20|) . we compute 

RiiL V) = Civ) cos{VSO + Div) siniVSO + Ripi^, v) , 
RipiC^ v) = ci cos{3VS^) + C2 sin(3\/^C) 



22 C3 cos([kj — VS]£,) + C4 cos([kj + v^]^) + 4 sin([Kj — vS]£,) + c^ sm{[Kj + vS]^) 

(21) 



j=i 



where j G {1 , 2} and 



a 



2,26 ^ ^' 32(5 ^ 

^ 2kj{kj-2V6)' ^ 2kj(kj +2\/^) ' 



V3B ^ _ V,B 

2kj{kj - 2V6) ' ^ 2kj{kj + 2VS) 



4= . /'\r.. ^ 4- . /^'"„^. - (22) 



Inserting (|15|l and H21|) into Eq. (|20|l and expanding the resulting equation 
trigonometrically yields 19 harmonics (that are also present for sines), which 
we list in Tabic Q We indicate which of these harmonics are always secular, 
sometimes secular, and never secular. 

At this order of perturbation theory, one finds 2: 1 (primary subharmonic), 
4: 1 (secondary subharmonic), 1 : 1 (harmonic), 2:1-1-1 (additive ultrasubhar- 
monic), and 2:1 — 1 (subtractive ultrasubharmonic) resonances. The first three 
types of resonances can occur with respect to either ki or K2, whereas the latter 
two require the interaction of both superlattice wavenumbers. Harmonic and ul- 
trasubharmonic spatial resonances have not been analyzed previously for BECs, 
and subharmonic resonances have only been analyzed in the case of regular OL 
potentials. At 0(e), we considered the case without 2: 1 resonances, so the as- 
sociated resonance conditions {kj — ±2v5) are necessarily not satisfied at the 
present [0(e^)] stage, as indicated in Tabled Second-order subharmonic (4:1) 
resonances have been studied in BECs in regular OL potentials 021 HOI ■ Their 
associated resonance conditions are kj — ±4v(5. (We return to subharmonic 



Label 


Harmonic 


Secular? 


Resonance when secular 


1 


cos{VSS,) 


Yes 


N/A 


2 


cos{3VS0 


No 


N/A 


3 


cos(5V^0 


No 


N/A 


4 


cos([ki - VS]S,) 


Assumed not in resonance at 0(e) 


2:1 


5 


cos([ki + VS]C) 


Assumed not in resonance at 0(e) 


2:1 


6 


cos([k2 - VS]£,) 


Assumed not in resonance at 0(e) 


2:1 


7 


cos([k2 + VS]£,) 


Assumed not in resonance at 0(e) 


2:1 


8 


cos([ki - 3VS]£.) 


Sometimes 


4:1 


9 


cos([ki + SVS]^) 


Sometimes 


4:1 


10 


cos([k2 - 3VS]0 


Sometimes 


4:1 


11 


cos([k2 + 3V^]0 


Sometimes 


4:1 


12 


cos([2ki - VS]0 


Sometimes 


1:1 


13 


cos([2ki + \/^]e) 


Sometimes 


1:1 


14 


cos([2k2 - VS]0 


Sometimes 


1:1 


15 


cos([2k2 + VS]£,) 


Sometimes 


1:1 


16 


cos([«;i + K2 - \/^]0 


Sometimes 


2:1 + 1 


17 


COs([ki + K2 + vS]^) 


Sometimes 


2:1 + 1 


18 


COs([/Ci — K2 — \/^]0 


Sometimes 


2:1-1 


19 


COs([ki — K2 + \/^]0 


Sometimes 


2:1-1 



Table 1: The harmonics in the right-hand-side of Eq. (|20|l after the formulas for 
Ro H15|) and i?i H21|l are inserted. We only list the cosines in this table, but the 
sines of these harmonics are present as well. We designate which harmonics are 
always secular, sometimes secular (under an appropriate resonance condition, 
as detailed in the text), and never secular. 

resonances in the case of OSLs later when we apply Hamiltonian perturbation 
theory.) The resonance relations for harmonic resonances are Kj = zt^/^. We 
will consider solutions that have harmonic resonance with respect to the primary 
lattice wavenumber ki. The resonance relation for additive ultrasubharmonic 
resonances is K2 + '^i = ±2\/J, and that for subtractive ultrasubharmonic res- 
onances is K2 — Ki — ±2v<^. In the remainder of this section, we consider the 
non-resonant, harmonically resonant, and the two types of ultrasubharmonic 
resonant states in turn. 

It is also important to remark that with the slow spatial variable t] ~ ex, the 
approximate solutions R{x) obtained perturbatively are valid for \x\ < 0(e^^) 
despite the fact that we employ a second-order multiple scale expansion. By 
incorporating a third ("super slow") scale e'^x, which is more technically de- 
manding, one can obtain approximate solutions that are valid for |a;| < 0(e~^) 

M 

Before proceeding, we also remark that in light of KAM theory, one expects 
different dynamical behavior (at least mathematically) depending on whether 
K2/K1 is an integer, a rational number, or an irrational number. Only the 



10 






situation K2 = 3ki has been prepared experimentally, so we concentrate on that 
case in our numerical simulations. 

We note additionally that we simulated the dynamics and examined the 
stability of MAWs using a numerical domain with periodic boundary conditions. 
This allows us to handle integer or rational values of AC2/K1 with appropriate 
selection of the domain parameters (so that the box size is an integer multiple 
of both spatial periods). However, quasiperiodic potentials cannot be tackled 
numerically within this framework for the extended wave solutions considered 
in this section. Our analytical work on MAWs is valid for all real ratios K2/K1. 

3.1 The Non- Resonant Case 

In the non-resonant case, effective equations governing the 0(£^) slow evolution 
are 

^' ^ TTx ^ [(/i("; 5^ Ki,K2)B'^ + /2(a, S, Ki,K2)A'^ + fs^a, S, ki, K2, 61)) D + /^{a, S, ki,K2)ABC 

+h{a, S, Ki,K2)B^ + fe{a, S, ki,K2)A'^B^ + fjla, 5, ki,K2)A^B + fs{a, S, ki, K2, 62)^] , 
1 

+h{a, (5, Ki,K2)A^ + fe{a, S, ki,k2)A^B^ + fr{a, S, ki, k2)AB^ + Js,{a, 5, ki, ^2)^] , 

(23) 

where 

A((5, Ki, K2) = 256(5^/2 (16(5^ - A5kI - A5kI + kIkJ) (24) 

and 

/i(a, S, Ki, K2) = 3/2(0?, (5, Ki, K2) , 
/2(a, 6, Ki, K2) = 96a6[16S^ - 4(5(k? + nj) + kIkJ] , 
fsia, S, Ki,K2,bi) = 2b6S%[-KlKl + AS{kJ + nl) - 165^] , 
h{a, 5, Ki, K2) = 2/2(0, (5, Ki, K2) , 
/5(a,(5, Ki,K2) = 15o^[-165^ + 4(5(ki + kI) - kJkI] , 
feia, S, Ki, K2) = 2/5(0, 5, Ki, K2) , 
fria, S, Ki, K2) = /5(a, S, ki, K2) , 

f8ia,6,Ki,K2,b2) = 64S[V^kI + V^k\ - ^biy^ + V^ + k?«;262) + \%b%2{.K{ + kI) - %\&%2\ ■ 

(25) 

In this case, the OSL does not contribute to O(e^) terms. 
Equilibrium solutions of 123|) satisfy 

_ (/li?^ + hA^ + h){hA^ + !^A^B^ + hAB^ + /sA) - [j^AB){hB^ + /eA^i?^ + ^,A^B + /gi?) 

/IA252 _ (y^52 + ^2^2 + J3)(J^^2 + ^2^2 + J3) 
(/lA^ + /2i?^ + /3)(/5g' + !&A^B^ + AA^jj ^ ^^g) , (/4Ai3)(/5A5 + /gA^i?^ + i,AB^ + /sA) 

(26) 
11 
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Figure 1: Evolution of the non-resonant spatially extended solution of Eq. p2|l 
with C and D in Eq. H21|l given by Eq. H23I) [see text for parameter details] for 
an OSL potential with V2 = 2Vi = 2 and K2 — 3ki — 3. The left panel shows 
the spatio-temporal evolution of |u(a:,t)p by means of a colored contour plot. 
The right panel shows spatial profiles of |up at four values of time {t = 100, 
200, 300, and 400). 

where one inserts an equilibrium value of A and B from Eq. H19|l . One then 
inserts equilibrium values of A, B, C, and D into Eqs. (|15|1 and H21I) to obtain 
the spatial profile R = Ro+ eRi -I- 0(£^) used as the initial wavefunction in the 
numerical simulations of the full GP given by Eq. (j3Jl . 

A typical example of the non-resonant case is shown in Fig. ^ with V2 = 
214 = 2 and K2 = 3ki — 12\/S — 3n/{2b), where b is the stretching factor given 
by Eq. (|ll|l . In this simulation, we used bi ~ b2 ~ 1 and e = 0.1. It can be 
clearly seen that the relevant solution is dynamically stable, which we found to 
be generic in our numerical experiments. Simulations with rational K2/K1 reveal 
similar phenomena. 
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3.2 Resonances 

In this subsection, we consider harmonic resonances, additive ultrasubharmonic 
resonances, and subtractive ultrasubharmonic resonances. In the evolution 
equations for the slow dynamics, one inserts the appropriate resonance relation 
into A and /1-/7. The function /§ has both the non-resonant contributions 
discussed above and additional resonant terms due to the OSL. Note addition- 
ally that there is a symmetry-breaking in the resulting equations because the 
functional form of the lattice contains only cosine terms. 

3.2.1 Harmonic Resonances 

When Kj — ±vo, there is a harmonic resonance. The effective equations gov- 
erning the O(e^) slow evolution in the presence of a harmonic resonance with 
respect to the primary lattice wave number ki are 

C" = T-7 ;^ [(fi{a,Ki,K2)B^ + /2(a,Ki,K2)-4^ + /si", ^i, ^2, 61)) D + /^{a, ki, K2)ABC 

A(ki, k2) 



+h (a, '^i , K2)-B^ + /e (a, ki , ^2)^^^^ + hia, ^i , H2)A'^B + f^s (a, ki , K2 , fe2)-S] , 

-r? ?r: [(/i("'«i''*2)^^ + h{a,Ki,K2)B'^ + h{a, Ki, K2,bi)) C + /4(a, ki, K2 

A(ki, k2) 



where 
and 



(27) 
A(ki, K2) = 768k?(4k? - kI) (28) 

/l(a,Ki,K2) = 3/2(0;, Ki,K2), 

/2(a, Ki, K2) — 288aK^(K2 ~ 4k^) , 
h{a,Ki,H2,bi) = 768ki5i(-K2 + 4ki) , 

fi{a, Ki,K2) = 2/2(0, Ki,K2) , 

/5(q!,ki,K2) = 45a^(-K2 + 4ki), 

/6(a,Kl,K2) = 2/5(0,(5, Kl,K2), 
/7(o,Ki,K2) = /5(o,(S, Ki,K2), 
f8s{a,Ki,K2,b2) = /„o„(o,Ki,K2) -t-32Vf (^2 - 4k^) , 
/8c(o,Ki,K2) = /non(o,Kl,K2) " 160Vf (Kj - 4Ki) , 

/„o„(o, «i, K2) = 192K2(y| - 4KlKlb2 + 16^%) ■ (29) 

If considering a harmonic resonance with respect to the secondary lattice wavenum- 
ber K2, one obtains the appropriate equations for the 0{e^) slow evolution by 
switching the roles of ki and K2. Note that the form of equations H29I) corre- 
sponds to H25|l except for the extra terms in /gc and fss that arise from the 
superlattice. 
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Figure 2: Same as Fig. Q] but for the harmonic resonant case with respect to 
the primary lattice wavelength. The solution given by Eq. p2l) is used as an 
initial condition with C and D in Eq. H21|) given by Eq. H27() with the functions 
H28U29(I [see text for parameter details]. 

The equilibria of l|27|l are given by Eq. I|26|l except that one inserts the 
functions from H29|l . Additionally, the expressions for C and D have /ss rather 
than /g as a prefactor for B and fsc rather than /g as a prefactor for A. One 
also inserts an equilibrium value of A and B from Eq. p9() . One then inserts 
equilibrium values of A, B, C, and D into Eqs. H15() and H21() to obtain the 
spatial profile R — Rq + ei?i + O(e^) to use as an initial condition in direct 
numerical simulations of Eq. ||3Jl. 

A typical example of the single-wavelength resonant case is shown in Fig. 
121 with V2 = 2Vi = 2 and K2 — 4ki = 4^/5 = n/b, where b is the stretching 
factor of Eq. Hll() : we used 61 = 62 = 1 and e = 0.1. The resulting (spatial) 
quasiperiodic patterns were generically found to persist in the dynamics of the 
system as stable (temporally oscillating) solutions. 
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3.2.2 Ultrasubharmonic Resonances 

Studying BECs in an OSL rather than in a regular OL allows one to examine 
the ultrasubharmonic spatial resonances resulting from interactions between the 
two lattice wavelengths ^J. As with harmonic resonances, an 0{e^) calculation 
is required to perform the analysis. 

When K2 + Ki — ±2vJ, one has an additive ultrasubharmonic resonance. 
The effective equations governing the O(e^) slow evolution in this case are l|?7|l 
with 

A(ki, K2) = 32kiK2(k1 + 2k2)(2ki + K2)(k1 + K2f (30) 

and 



+ K2) + 34kiK2(«i + 1^2) + 46K1K2] 
/5(a, Ki, K2) = 15a^KiK2[5KiK2 + 2{kI + nl)] , 

/6(a, Ki,K2) = 2f5{a,S,Ki,K2), 
/?(", Ki,K2) = f5{a,S,Ki,K2), 
f8sia,Ki,K2,b2) = fnon{a,Ki,K2) - /res («,«!, K2) , 
fsc{a,Ki,K2,b2) = fnon{a,Ki,K2) + fres{a, Ki, K2) , 

fnonia, Ki, K2) = 16[13k^K262(Ki + Kj) + 46KiK2fo2 + Sk^KjI^I^ + ^2^) 

+ 2KlK2(Vf K? + V^KJ + Kfb2 + K^) + 34k?K^62(k? + K^) 

+ AK^K2{V^nl + Vi^l) + V^4 + Vi4] , 
/res(a, Ki, K2) = 32ViV2[7k?k2 + (k^ + 4) + 4kiK2(k? + k^)] • (31) 

Note that all the terms in fres are proportional to VlV2, as they arise from the 
effects of interacting lattice wavelengths. 

Equilibria in this situation again satisfy 1)26(1 except that one now inserts 
functions from H30II31() . Again, the expressions for C and D have /ss rather 
than /g as a prefactor for B and /gc rather than /§ as a prefactor for A. One 
again inserts an equilibrium value of A and B from Eq. ((19|l . One then inserts 
equilibrium values of ^, B, C, and £> into Eqs. ((15|l and H21|l to obtain the 
initial spatial profile i? = i?o + e^i + O(e^). 

A typical simulation of an ultrasubharmonic resonance is shown in Fig. 13 
with V2 = 2Vi = 2 and K2 = 3ki = 3vo/2 = 37r/(86), where b is again given 
by Eq. ((11(1 with hi — b2 — 1 and e = 0.1. The resulting complex patterns were 
found to persist as stable dynamical structures (with periodic time dynamics). 

When K2 ^ '^1 = ±2v<^, one has a subtractive ultrasubharmonic resonance. 
The effective equations governing the 0{e^) slow evolution in this case are again 
(EH), with 

A(ki, K2) = 32kiK2(ki — 2k2)(2ki — K2)(ki — K2)^ (32) 
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Figure 3: Same as Fig. ^ but for an additive ultrasubharmonic resonance, 
which arises from the interaction of the EEC's two wavelengths. The solution 
of Eq. (|12|l is used as an initial condition with C and D in Eq. H21|l given by 
Eq. (|27|l with the functions 1)301131(1 [see text for parameter details]. 
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and 

/l(a,Ki,K2) = 3/2(a,Ki,K2), 

/2(a, Ki, K2) = 24aKiK2[-2(Kf + K2) + 9(k? + k^) ~ Unjnj] , 
fsia, Ki, K2, bi) = 16kiK2&i[2(ki + K2) - l3KiK2{Kt + ^2) + 34:kIkI{kI + nl) - 46KiKf] 
/4(a, Ki, K2) = 2/2(a, Ki, K2) , 
/5(a, Ki, K2) = 15a^KiK2[-5KiK2 + 2(k^ + K2)] , 
/6(a, Ki, K2) = 2/5(0;, (5, Ki, K2) , 

/7(a, Ki,K2) == /5(a, (5, Ki,K2), 
/8s(a, Kl,K2,62) = fnonia,Ki,K2) - fres{a, Ki, K2) , 
/8c(a, Kl,K2,^2) = fnon{a,Ki,K2) + /res (",«!, ^2) , 

fnonia, Kl, K2) = 16[-13KiK2^2(Ki + ^2) - 46KiK2^2 " ^1^X1^2(^1 + ^2^) 

+ 2kiK2(V'2^K? + V"i2«-| + K?fc2 + ^^62) + 34KfK|62(K? + K^) 

+ 4Ki«2(n'«? + V^nl) - V^4 - Vi^] , 

fres{a,Ki,K2) = 32^1^2 [-7^1 Kj - (Kj + K2) +4kiK2('«i + Kj)] • (33) 

As with the additive uhrasubharmonic resonance, all the terms in f^es are pro- 
portional to VlV2. 

Equilibria in this case again satisfy H26|l except that one inserts the functions 
from (|32II33|I . Recall once more that the expressions for C and D have fss rather 
than /g as a prefactor for B and fsc rather than /g as a prefactor for A. One 
also inserts an equilibrium value of A and B from Eq. (|19|l . One then inserts 
equilibrium values of A, B, C, and D into Eqs. p5|) and (|21|l to obtain a spatial 
profile R = Rq + eRi + O(e^) to utilize as an initial wavefunction in numerical 
simulations of Q. In this case, the numerical simulations yielded similar (stable) 
temporal dynamics as for additive ultrasubharmonic resonances. 

4 Hamiltonian Perturbation Theory and Sub- 
harmonic Resonances 

In this section, we build on recent work |49l I5(JJ and apply Hamiltonian pertur- 
bation theory to pi)|l to examine period-multiplied wavefunctions and spatial 
subharmonic resonances in repulsive BECs loaded into OSL potentials. (For 
expository reasons, we repeat some details of the derivation from those works 
in the present one.) We perturb off elliptic function solutions of the underly- 
ing integrable system and study 2n : 1 spatial resonances with a leading-order 
perturbation method. Perturbing off simple harmonic functions, by contrast, 
requires a perturbative method of order n to study 2n : 1 resonances. At the 
center of KAM islands lie 'period-multiplied' states. When n — 1, one obtains 
period-doubled states in u corresponding to 2 : 1 subharmonic resonances. Our 
analysis reveals period-multiplied solutions of the GP (Q with respect to both 
the primary and secondary sublattice. 
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The dynamical systems perspective on period-doubled states and their gen- 
eralizations for BECs in OSL potentials given here complements theoretical 
and experimental work by other authors for the case of regular OL potentials. 
In recent experiments, Gemelke et al. ,26] constructed period-doubled wave- 
functions, which have received increased attention (for regular lattices) during 
the past two years. In earlier work, Smerzi et al. |56j reported theoretical 
studies of spatial period-doubling in the context of modulational ( "dynamical" ) 
instabilities of Bloch states and Cataliotti et al. jj^l reported experimental ob- 
servations of superfluid current disruption in chains of weakly coupled BECs. 
Period-doubled states, interpreted as soliton trains, then arise from dynamical 
instabilities of the energy bands associated with Bloch states [321 • 

4.1 Unforced Duffing Oscillator 

We employ exact elliptic function solutions of Buffing's equation [Eq. H10|l with 
Vi = V2 = 0] , so we no longer need to assume the coefficient of the nonlinearity 
is small. Therefore, we use the ODE 

R" + SR + aR^ + eR[Vi cos{kix) + V2 cos(K2a;)] = , (34) 

which is just like (|l()|l except that a no longer has the prefactor e. 

When e = 0, solutions of ()34|l are expressed exactly in terms of elliptic 
functions (see, e.g., |H21|Sni and references therein): 

R = apcn{u,k) , (35) 



where 



uix + uq , Ui = S + ap . 



e "^' 



2(<5-f ap2) ' 
ui > 0, p>0, k^e^, ere {-1,1}, (36) 

and mq is obtained from an initial condition (and can be set to without loss 
of generality). When ui G M, the solutions given by l|5S)l are periodic. When 
A;2 < 0, which is the case for repulsive BECs with positive chemical poten- 
tials, equation H36|) is interpreted using the reciprocal complementary modulus 
transformation (as discussed in Ref. (5U|V 

Equation (|34|l is integrated when e = to yield the Hamiltonian 



with given energy 



2 2 4 



i'''(^^ + '"''>4(T^' P«' 
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where fc'^ := 1 — k^. 

The center at (0, 0) satisfies h — p'^ — P — 0. The saddles at {±^/—6/a, 0) 
and their adjoining separatrix (consisting of two heterochnic orbits) satisfy 

h^-^, p^ = ±, k^^-oo. (39) 

The sign a — +1 is used for the right saddle and cr = — 1 is used for the left one. 
Within the separatrix, all orbits are periodic and the value of a is immaterial. 

4.2 Action-Angle Variable Description and Transforma- 
tions 

For the sake of exposition, we construct an action-angle description in steps. 
First, we rescale (|34|) using the coordinate transformation 

X^VSx, r=.L^R (40) 

V 

to obtain 

r" + r - r^ = (41) 

when Vi = V2 = 0. In terms of the original coordinates, 

a 



R{x) = \ r(V5x\ . (42) 



The Hamiltonian corresponding to (|41|) is 

Ho{r,s)^^s^ + ^r'-^r^^h, /i € [0, 1/4] , (43) 

where s := r' = dr/dx- Additionally, p^ G [0, 1] and 

With the initial condition r(0) — p, s(0) = (which implies that uq = 0), 
solutions to H41() are given by 

six)^-p[l-p'f'^^{[^-p']'^'x,k)dn[[l~p^]'/\,k) . (45) 

The period of a given periodic orbit F is 

T{k) =idx= ^fflf , (46) 
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where 4:K{k) is the period in u of cn(u, k) [SH]- The frequency of this orbit is 



n{k)^ 



2K{k) 



(47) 



Let Th denote the periodic orbit with energy h = Ho{r,s). The area of 
phase space enclosed by this orbit is constant with respect to x, so we define 
the action ^7\ 

If 1 r'^C') 



'■■=2. 
which is evaluated to obtain 



sdr = 



2tt 



5(x)] dx, 



J=^^^^[E{k)-{l-py2)Kik)] 
The associated angle in the canonical transformation (r, s) 



(J,$)is 



(48) 



(49) 



(50) 



The frequency il.{k) decreases monotonically as k^ goes from — cx3 to [that 
is, as one goes from the separatrix to the center at {r,s) = (0,0)]. With this 
transformation, equation 1451) becomes 

r{J, $) = p{J) en (2if (fc)$/7r, k) , 

six) = -p( J) Vl - PiJ)^ sn {2K{k)<^/7r, k) dn (2i^(fc)$/7r, k) , (51) 

where k = k{J). 

After rescaling, the equations of motion for the forced system (|34|l take the 
form 

„3 



r + r — r 



y,cos(^X^ 



y.cosl^^X 



r = 



(52) 



with the corresponding Hamiltonian 
H{r,s,x) = Ho{r,s) +eHi(r,s,x) 

In r Q ^ A £ 



2 I -'^ 2 -'^4 

-s + —r -r 
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F,cos(^X^ 



1^2 cos (^X^ 



(53) 



In action-angle coordinates, this becomes 

iJ($, J, x) = \p{Jf - \p{Jf + ^P( J)' cn^ (2if (fc)$/7r, k) 



l/,cos(^x)+T/.cos(^X 

(54) 



A more convenient action-angle pair (^,j) is obtained using the canonical 
transformation ($, J) — > ('/', j), defined by the relations 



j{J)^\p{J)\ ^^^^^) = J(f) 



(55) 
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where 



e = 



2j-l' 



J{j) = 37^^27 [e{j) - (1 - j)K{j) 
K{j) = -K[k{j)] , E{j) = -E[k{j)] . 



Additionally, 



dJ 



l-2i 



Aj)^-f-VT^KU) ^^^^ 



(56) 



(57) 



Note that J ^ j for small-amplitude motion. Furthermore, j = at the origin, 
and J = 1/2 on the separatrix. The Hamiltonian H54|l becomes 



-ff (0, Ji X) = J - f + TJ cn^ 



Kill 
J'U)' 



Vi cos f -j=x ] + V2 cos 



K2 



(58) 



4.3 Perturbative Analysis 



A subsequent 0{e) analysis at this stage allows us to study 2n: 1 subharmonic 
resonances for all n G Z. Fourier expanding the en function yields 



en 



Kjj) 



,k\^Bo{j)+f2Bicos(^^^ , (59) 



where the coefficients Bi{j) are obtained by convolving the Fourier coefficients 

iHaEm, 



BnU) 



k{J)K{j) 



bn[k{j)] 



b„{k) = -sech [{n + 1/2) TrK'{k)/K{k)] , 



(60) 



of the en function in Eq. (|58() . where K'{k) := K{\^1 — fc^) is the complementary 
complete elliptic integral of the first kind j59l ^ . 

The resulting 0{e) term in the Hamiltonian (|58() is 



£Hi{(f>,j,x) ^ -^jBoiJ) 



F,cos(^x)+^.cos(^X 






2l(j) Ki \ / 2l(j) m 

cos I ^TT-r H -X + cos — -7— ^x 



J'U) V6' 



J'U) V6' 



l'=l 



2l'q 



J'U) ^fs' 



K2 

—j=X I + COS 



2l'a 



K2 



J'U) Vs' 



(61) 
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The Hamiltonian (|61|l is an expansion over infinitely many subharmonic 
resonance bands for each of the primary and secondary lattice wavenumbers. 
Each resonance corresponds to a single harmonic in Ht)l|) . To isolate individual 
resonances, we apply the canonical, near-identity transformation |62[ I5UJ 

to H61() with an appropriate generating function Wi that removes all the res- 
onances except the one of interest. The subscript i in Qi designates whether 
one is considering a resonance with respect to the primary or secondary lat- 
tice wavenumber. The transformation Ht)2|l is valid in a neighborhood of this 
2n : 1 resonance and yields an autonomous 1 DOF resonance Hamiltonian that 
determines its local dynamics, 

K{Q, P,x;n)^P-P' + ^V.PBr^P) cos (^|||^ - ^x) + 0{e^) . (63) 

In focusing on a single resonance band in phase space, one restricts P to a 
neighborhood of P„ , which denotes the location of the nth resonant torus asso- 
ciated with periodic orbits in 2n: 1 spatial resonance with the primary (i = 1) 
or secondary {i = 2) sublattice (recall that ki < ^2). 

The resonance relation associated with 2n : 1 resonances with respect to the 
ith sublattice is (SUI 

^ = ±2nn(P„) . (64) 

vo 
Because J7 < 1 is a decreasing function of P e [0, 1/2), the associated resonance 
band is present when 

^ < 2n . (65) 

For example, when Ki — 2.5 and <^ = 1, there are resonances of order 4:1, 6:1, 
8:1, etc, but there are no resonances or order 2:1. Analytical expressions for 
the sizes of the resonance bands and the locations of their saddles and centers 
are the same as those obtained for BECs loaded into OLs; they are derived in 

Ref. inni- 

To examine the time-evolution of period-multiplied solutions, we need only 
the locations of centers, which are obtained by applying one more canonical 
transformation. We use the generating function 

F,iQ,,Y,x;n) = Q,Y-^^JiY)x, (66) 



which yields 



oQi 

r)P K 

i = ff^iQ^^y^x) = Q^- ^J'{Y)x. (67) 
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The resonance Hamiltonian (|63|l becomes 

dF 

K^i^,Y) ^ K{Q,,P,X;n) + -^iQ,,Y,x) 

ox 

which is integrable in the [Y, S,) coordinate system. 

The centers of the KAM islands associated with this resonance occur at [HO] 

Y, = r„ + eAY + 0{e^) , (69) 

where 

1 r R (V \ Ji^v R' (V \ ^ 

(70) 



^25 



Bn{Yn)+YnB'jY^) 



_r!(y„)vi-2r„i^'(r„)-i. 

and the sign is — when n is even and + when n is odd. One then converts the 
value Yc back to the original coordinates to obtain an estimate [Rc^ Sc) of the 
location of the center in phase space. [One obtains the locations of the other 
centers associated with the same resonance band using iterates of {Re, Sc) under 
a Poincare map, but we only need one of these centers for a given resonance to 
examine the time-evolution under the GP equation ^ of these solutions, which 
provide the initial wavefunctions for the PDE simulations.] 

In our numerical computations, we use the parameter values h — 2m — 1, 
S = I, a = -1, e = 0.01, and Vl = 1 in Eqs. POI). With k = 1.5, there is a 
center for the 2 : 1 resonance with respect to the primary sublattice at Re ~ 0.753 
and Sc = 0, so one uses R = 0.753 cos(Kia:/2) as an initial wavefunction in simu- 
lations of © for any height V2 and wavenumber K2 of the secondary sublattice. 
Such a solution is shown in Fig. 0] for V2 = 2 and K2 = 3ki. It is dynamically 
stable and sustains only small amplitude variations (but is otherwise essentially 
stationary). One can similarly examine initial wavefunctions corresponding to 
2 : 1 resonances with respect to the secondary sublattice. 

With Ki = 2.5, there is a center for the 4 : 1 resonance with respect to 
the primary sublattice at {Re, Sc) ~ (0.691,0.324), so (recalling the chain rule) 
one uses R — 0.691 cos(Kia;/4) -I- 0.518 sin(Kix/4) as an initial wavefunction 
in simulations of (O The results with K2 = 3ki and V2 = 2 are shown in 
Fig. We observe a wriggling pattern in the contour plot (in the left panel), 
which indicates (structurally stable) spatio-temporally oscillatory behavior of 
the condensate. 

With Ki = 3.8, there is a center for the 6 : 1 resonance with respect to the 
primary sublattice at Re ~ 0.859 and Sc — 0, so one uses R = 0.859 cos(kix/6) 
as an initial wavefunction in simulations of ((2Jl. We observe that this period- 
multiplied state is stable with small- amplitude oscillations, as was the case for 
2 : 1 resonances. At the same value of ki, there is a center for the 8 : 1 resonance 
with respect to the primary sublattice at Re ~ 0.9354 and Sc ~ 0.0718, so one 
uses R — 0.9354 cos(Kia;/8)-t-0.151 sin(Kia;/8) as an initial wavefunction. As was 
the case for 4 : 1 resonances, PDE simulations reveal structurally stable spatio- 
temporally oscillatory behavior of the condensate (shown for 4 : 1 resonances as 
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Figure 4: Same as Fig. ^ but for a 2 : 1 resonance with respect to the primary 
sublattice. The solution described in the text [R = 0.753 cos(Kia:/2) with ki = 
1.5 = K2/3] is used as the initial condition (see the text for further parameter 
details). The solution appears to be dynamically stable and only sustains a 
small-amplitude oscillation. 
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Figure 5: Same as Fig. Q] but for a 4 : 1 resonance with respect to the pri- 
mary sublattice. The solution described in the text [R — 0.691 cos(Kia;/4) -|- 
0.518 sin(Kia;/4) with ki = 2.5 = K2/3] is used as the initial condition (see 
the text for further parameter details). While structurally stable, the solution 
pattern appears to be a wriggling one, indicating a spatio-temporal breathing. 
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a wriggling pattern in the left panel of Fig.EJ. This difference between "odd" 
and "even" subharmonic resonances arises from the fact that the former contain 
centers on the i?-axis, whereas the latter do not. The resulting initial conditions 
in the even case hence require both sine and cosine harmonics, resulting in the 
observed spatio-temporal breathing. 

From a more general standpoint, resonance bands emerge from resonant 
KAM tori at action values P* that satisfy a (three-term) resonance relation 
with respect to both sublattices |CTICT) . 

m^+n2-^ = 2nr!(n), (71) 

where n, rii, and 712 all take integer values. The single-sublattice resonance 
relation l|M|l is a special case of (|7T]| . 

5 Conclusions 

In this work, wc analyzed spatially extended coherent structure solutions of the 
Gross-Pitaevskii (GP) equation in optical superlattices describing the dynam- 
ics of cigar-shaped Bose-Einstein condensates (BECs) in such potentials. To 
do this, we derived amplitude equations governing the evolution of spatially 
modulated states of the BEG. We used second-order multiple scale perturba- 
tion theory to study spatial harmonic resonances with respect to a single lattice 
wavenumber, as well as additive and subtractive ultrasubharmonic resonances. 
Harmonic resonances are a second-order effect that can occur in regular peri- 
odic lattices, but ultrasubharmonic resonances can only occur in superlattice 
potentials, as they arise from the interaction of multiple lattice wavelengths. In 
each situation, we determined the resulting dynamical equilibria, which repre- 
sent spatially periodic solutions, and examined the stability of the corresponding 
solutions via direct simulations of the GP equation. In every case considered, 
the solutions (non- resonant, resonant with a single wavelength, and resonant 
due to interactions between two wavelengths) were found numerically to be dy- 
namically stable under time-evolution of the GP equation. Finally, we used 
Hamiltonian perturbation theory to construct subharmonically resonant solu- 
tions, whose spatio-temporal dynamics we illustrated numerically in a number 
of prototypical cases. 
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